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Planets opening dust gaps in gas disks 
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Abstract. We investigate the interaction of gas and dust in a protoplanetary disk in the presence of a massive planet using a 
new two-fluid hydrodynamics code. In view of future observations of planet-forming disks we focus on the condition for gap 
formation in the dust fluid. While only planets more massive than 1 Jupiter mass (Mj) open up a gap in the gas disk, we find 
that a planet of 0. 1 Mj already creates a gap in the dust disk. This makes it easier to find lower-mass planets orbiting in their 
protoplanetary disk if there is a significant population of mm-sized particles. 
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1. Introduction 

Resolved images of protoplanetary disks can reveal density 
structures in gas and dust hinting at the presence of embed- 
ded planets. It is therefore important to investigate the re- 
sponse of a dusty circumstellar disk on planets, as opposed to 
for example radiation pressure effects in optically thin disks 
jTakeuchi & Artvmowiczl boo ill . 

The most dra matic disk structure a planet can create is an 
annular gap ( Papal oizou & Linll984l) . and these planet-induced 
gaps can be used as t racers for recent giant planet formation. 
Recently, IWolf et all (j2002) showed that the Atacama Large 
Millimeter Array (ALMA) will be able to detect a gap at 5.2 
AU caused by the presence of a Jupiter mass planet in a disk 
140 pc away in the Taurus star-forming region. However, they 
assumed a constant dust-to-gas mass ratio of 1:100, whereas in 
the presence of large pressure gradien ts the gas and the dus t 
evolution will decouple to some extent (Takeuchi & Lin 2002). 
Due to radial pressure support the gas orbits at sub-Keplerian 
velocity, and when the dust particles are slowed down by drag 
forces they move inward. 

In this letter we present the first results of two-dimensional 
numerical simulations where we treat the gas and the dust as 
separate but coupled fluids. These two-fluid calculations take 
into account the full interaction of gas and dust, and can be 
used to simulate observations in a better way. 

In Sect. |2 we briefly discuss the physics of gas-dust inter- 
action. We describe the numerical method in Sect. |5] and the 
initial conditions in Sect. 0] We present the results in Sect. |5] 
and we conclude in Sect.|6] 



2. Basic equations 

Protoplanetary disks are fairly thin, i.e. the vertical thickness 
H is small compared with the distance r from the centre of 
the disk. Typically we use H/r — h = 0.05. Therefore it is 
convenient to average the equations of motion vertically, and 
work with vertically averaged quantities only. The governing 
equations as well as the gas disk model are the ones described 
in detail in lKlevNl999l) . 

2.1. Dust 

We treat the dust as a pressureless fluid, and its evolution is 
governed by conservation of mass and conservation of lin- 
ear and angular momentum. The major difference between the 
equations for the gas and the dust is the absence of pressure in 
the latter. In this sense the dust fluid behaves like a gas that is 
moving always with supersonic velocity. This implies that near 
shock waves, where the gas goes from sonic to supersonic flow 
and large density and velocity gradients are present, dust and 
gas behave in a very different way. 

The interaction between the gas and dust fluids occurs only 
through drag forces. The nature of the drag force depends on 
the size of the particles with respect to the mean free path of the 
gas molecules. We consider only spherical particles with radius 
s, and for the particle sizes we are i nterested in (s » 1 mm) 
we ca n safely use Epsteins drag law ( Take uchi & Artvmowiczl 
2001). In the limit of small relative velocities of gas and dust 
compared to the local sound speed, the drag forces are written 
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where Q,k is the Keplerian angul ar velocity, and T R is the non- 
dimensional stopping time ( iTakeuchi & Linl2 002): 
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Here, p p is the particle internal density, s is the radius of the 
particle and c s is the sound speed. We have used p p — 1.25 
g cm -3 , and the gas density p g is found by using p g = E g /2if. 
During simulations we have always checked whether Eq.^ap- 
plied. 

3. Numerical method 

For the evolu tion of the gas component, we use the method 
described by Paar dekooper & Mellemal J2004I) . It is a second- 
order Eulerian hydrod ynamics co de which uses an approxi- 
mate Riemann solver (Roe 198 ! ])■ It derives fr om a general 
relativistic method ( Eulderin k & Mellemalfl99l) . from which 
it inherits the capability to deal with non-Cartesian, non- 
inertial coordinate systems. We work in a cylindrical coordi- 
nate frame (r, (j>), corotating with the embedded planet that 
has a Keplerian angular velocity. A module for Adaptive Mesh 
Refinement (AMR) can be used to obtain very high resolu- 
tion near the planet. We use basically the same method to 
follow the evolution of the dust, but without the pressure 
terms. This is a very di fferent approach than for example in 
iKlahr & Henningl ill 9971) . who solve the equations of motion 
of particles instead of the fluid equations. All the source terms 
except t he drag forces are integrated using stationary extrap- 
olation ( Paardekooper & Mellemal l2004l) . The drag forces are 
incorporated separately using an ordinary differential equation: 



(3) 



Separating the drag forces from the advection in this way 
puts constraints on the time step. Basically the assumption 
is that no drag force acts on the dust during the advection, 
and this is clearly not correct when the drag force source 
terms are very stiff. The time step is set by the hydrodynamic 
Courant-Friedrichs-Lewy (CFL) condition, and therefore one 
has to make sure that the typical time scale for the drag forces 
(the stopping time) is not much smaller than this time step. 
Switching the order of integration every t ime step yield s a sec- 
ond order accurate update for the state JStrandfl968l) . If the 
stopping time is very small (well-coupled particles) the dust 
can not be assumed to move free of drag during a dust update. 
This puts a lower limit on the particle size that we can put into 
the simulations. For a typical disk, this limit is about 0.1 mm. 
The time step can however be reduced to include smaller parti- 
cles. 



4. Initial and boundary conditions 

The standard numerical resolution is (n r ,«0) = (128,384). 
This way, the cells close to the planet have equal size in the 
radial and in the azimuthal direction. We do not resolve the 
Roche lobe at this resolution, but since we are interested only 
in the large scale evolution of the disk, this is not important. 



We will discuss the flow within the Roche lobe and accretion 
processes in future papers. Using a low resolution allows us 
to run longer simulations, which is important because the drift 
velocities can be very low. We take a constant initial surface 
density So = 34 g cm~ 2 . This surface density corresponds 
to the value at 5.2 AU in a disk with 0.01 M Q within 100 
AU and a surface density profile that falls off as r -1 ' 2 , or to 
about 13 AU in the Minimum Mass Solar Nebula. The initial 
dust-to-gas ratio is 0.01. Using the radius of the planet's or- 
bit as our unit distance, the inner boundary is at r — 0.4 and 
the outer boundary at r = 2.5. The gas rotates with a slightly 
sub-Keplerian velocity due to the radial pressure gradient. The 
dust rotates exactly with the Keplerian velocity, initially. All 
radial velocities are take n to be zero. The constant kinematic 
viscosity is set so that the lShakura & SunvaevMl973l) a param- 
eter equals 0.004 at r = 1. For a disk described above, the 
critical planet mass for gap opening is about 1 Jupiter mass 
(1 Mj) dBrvden et alJll999h . We varied the mass of the planet 
between 0.00 1 M.t and 0.5 Mj . The boundary conditi ons are 
as in [ Paar dekooper & Mellema (2004), deriving from iGodonl 
( 1996). They are non-reflecting, so that all waves generated in 
the simulated region leave the computational domain without 
further influencing the interior. 



5. Results 

The evolution of the gas compo nent is described in detail in 
IPaardekooper & Mellemal d2004l) . and the drag forces on the 
gas are too small to alter the flow. Therefore we focus on the 
evolution of the dust particles. 

5.1. General flow structure 

For the standard case we use a planet of 0.1 Mj, and we con- 
sider dust particles of 1 mm. In the left panel of Fig.^ we can 
see that the dust particles are opening a gap, which leads to 
variations in the dust-to-gas ratio of more than an order of mag- 
nitude (middle panel). For this small planet these variations in 
dust density are so large that the spiral waves are no longer vis- 
ible. The radial cut (right panel) shows that this planet is not 
able to open up a gap in the gas, while the dust reacts more 
strongly and shows a very dense ring at r = 1.3 just outside a 
very low density gap. 

Also shown in the right panel of Fig. [I] are positions of the 
3:2, the 2: 1 and the 3:1 mean motion resonances. Especially the 
3:2 resonance (at r = 1.25) plays a role in structuring the disk, 
while only faint features can be seen at the other resonances. 
There also exists a small density bump at corotation (r = 1.0). 

Figure|2]shows a close-up on the gas density near the planet 
in a high-resolution simulation. The high-density core can be 
seen, as well as the spiral shocks. The arrows indicate the rel- 
ative velocity of the dust compared to the gas, showing that 
the dust decouples from the gas near the spiral shocks. In the 
shocks the gas velocity is reduced, and the gas is also deflected. 
The dust particles do not feel the shock directly, but notice its 
presence by the drag forces in the post-shock region. Figure [2] 
shows that it takes a certain distance before the velocity differ- 
ence between dust and gas has disappeared again. Because of 




Fig. 1. Dust flow near a 0.1 Mj planet. Left panel: Grey-scale plot of the logarithm of the dust surface density after 400 planetary 
orbits. Middle panel: Logarithm of the dust-to-gas ratio. Right panel: radial cut at <f> = (opposite to the planet). Solid line: 
gas surface density, dashed line: dust surface density x 100. The filled circles indicate the 3:2, the 2:1 and the 3:1 mean motion 
resonances. 



this, the flux of dust particles streaming into the spiral wave is 
larger than in the case of a perfectly coupled gas-dust mixture. 
Since the wave is responsible for gradually pushing the gas and 
dust particles away from the orbit of the planet this leads to 
a depletion of dust particles, and hence a deeper gap in the 
dust distribution. The dust density is affected most where the 
waves are the strongest, which is approximately at a radial dis- 
tance h from the orbit of the planet. This is where a low density 
gap starts to form, leaving the particles at r = 1 largely unaf- 
fected. The particles moving to larger radii seem to get trapped 
in the 3:2 and the 2:1 mean motion resonances, creating the 
dense rings at r = 1.3 and r = 1.6. Interestingly, it has been 
found that the Kuiper Belt Objects also tend to accumulate at 
the the 3:2 and 2: 1 mean motion resonances with Neptune (Luu 
& Jewitt 2002). Whether this is related to the effect we find, 
requires further study. Inside corotation, the combined inward 
force due to the spiral waves, the imposed radial temperature 
gradient and the viscous gas flow is so large that basically all 
particles move over the resonances quickly and end up near the 
wave-damping boundary where they are absorbed completely, 
leaving an almost featureless inner disk. Simulations with dif- 
ferent planetary masses showed that the minimum planetary 
mass for this mechanism to operate is 0.05 Mj with particles 
of 1 mm. 

5.2. Simulated observations 

To investigate the observational effects of the gas-dust decou- 
pling we simulate an image taken at a wavelength of 1 mm. The 
disk is optically thin at this wavelength, and we put the planet at 
5.2 AU. Calculations with different particle sizes showed that 
the minimum grain size for gap formation is approximately 0.1 
mm, with a very sharp transition. We therefore assume that all 



particles larger than this critical size create a gap, and smaller 
particles follow the gas exactly. To estimate the relative amount 
of mass in the se two families of particles we take an MRN size 
distribution ( Mathi set alJI 19771) . in which the number density 
of particles is a power law in s with exponent —3.5. For a maxi- 
mum particle size of 5 mm most mass is in particles larger than 
our critical size, and as a conservative assumption we take an 
equal amount of mass in the particles smaller than 0.1 mm and 
in the larger particles. But since t he opacity at 1 mm is domi - 
nated by particles with s m 1 mm ( Mivake & Nakagawa 1993), 
the emission will be dominated by the lar ger particles. We used 
the siz e- dependent opacity data from Mivake & Naka gawal 
(1993) to calculate the flux densities. In Fig. [5] the resulting 
flux densities are shown, convolved with a Gaussian of FWHM 
2.5 AU, which corresponds to an angular resolution of 12 mas 
at 140 pc (the Taurus star forming region). This resolution is 
comparable to the maximum resolution of 10 mas that ALMA 
will achieve. In the left panel of Fig. [5] we computed the flux 
under the assumption that all particles mov e with the gas, as 
was done for example bv IWolf et alJ J2002I) . There is a little 
dip visible due to the decrease in density in the gas, but only 
when we include the density distribution of Fig^for the larger 
particles a clear gap emerges (middle panel of Fig.|5Jl. The right 
panel of Fig.[3]shows the radial dependence of the flux density. 
From this panel we can see that the contrast between the inside 
of the gap and the gap edges is enhanced from almost a factor 
of two (dashed line) to more than an order of magnitude (solid 
line). So the impact of the dynamics of the larger particles on 
the appearance of a low-mass planet in a disk is considerable. 
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Fig. 3. Logarithm of flux densities at 1 mm, normalized by the maximum and convolved with a Gaussian of FWHM 2.5 AU, 
corresponding to a resolution of 12 mas at 140 pc. Left panel: all particles follow the gas exactly (static dust evolution). Middle 
panel: particles larger than the critical size decouple from the gas (dynamic dust evolution). Right panel: the corresponding radial 
flux densities. 



6. Summary 

We have presented the first two-fluid calculations of a circum- 
stellar disk with an embedded planet. These models allow us to 
make predictions on the dust distributions that will be observed 
in protoplanetary disks in the near future. 

For a reasonable disk mass (0.01 M within 100 AU) the 
strong spiral shocks near the planet are able to decouple the 
larger particles (> 0.1 mm) from the gas. This leads to the for- 
mation of an annular gap in the dust, even if there is no gap in 
the gas density. Because the opacity at millimeter wavelengths 
is dominated by these larger particles, the signatures of low- 
mass planets in disks can be stronger than previously thought. 
The minimum mass for a planet to open a gap this way was 
found to be 0.05 Mj for 1 mm particles. 
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Fig. 2. Close-up on the gas density near the planet for a high 
resolution AMR simulation (3 levels of refinement, so a reso- 
lution that is 8 times higher than our standard resolution). The 
arrows indicate the relative velocity of the dust compared to the 
gas. The typical relative velocity across the shocks is 0.035 c s . 
Note that the color scale is reversed with respect to Fig.^ 



